/*
 * Solver.h
 *
 *  Created on: 21.08.2010
 *      Author: jagiella
 */

#ifndef SOLVER_H_
#define SOLVER_H_

#include "SparseMatrix.h"

#define COMPARE_WITH_EXACT_SOLUTION 0

template <class T> class Solver
{
	//
	T *v0d;
	T *v1d;
	T *v2d;
	T *v3d;
	T *v4d;
	T *v5d;
	T *v6d;

public:
	// Solver Type
	static const char CG 			= 1;
	static const char BiCGSTAB		= 2;
	static const char Jacobi		= 3;
	static const char GaussSeidel	= 4;
	static const char SOR			= 5;
	static const char ILU			= 5;

	// Solver Type (Preconditioner)
	static const char Diagonal			= 6;
	static const char Tridiagonal		= 7;
	static const char LowerTriangular	= 8;
	static const char UpperTriangular	= 9;

	char 	type;
	int  	maxit;
	T	 	maxerr;
	int		size;
	T		omega;

	// Exact Solution
	bool 	useExactSolution;
	T*		exactSolution;

	Solver( int matrix_size, char type);
	~Solver();
	int Solve( SparseMatrix<T> *A, T *b, T *x);
	//int Solve( SparseMatrix *A, double *b, double *x);
	//int Solve( SparseMatrix *A, float *b, float *x);

	void calculate_x( T *x, T omega, T *w, T alpha, T *p, int N);
	void calculate_r( T *r, T omega, T *v, T alpha, T *s, int N);
	void calculate_p( T *x, T *y, T alpha, T *a, T beta, T *b, int N);
	void calculate_p1( T *x, T *y, T alpha, T *a, int N);
	void calculateScaledSum( T *o, T s1, T *v1, T s2, T *v2, T s3, T *v3, int N);

	int SolveConjugateGradient( SparseMatrix<T> *A, T *b, T *x);
	int SolverPreconditionedConjugateGradient( SparseMatrix<T> *A, SparseMatrix<T> *M, T *b, T *x);
	int SolveBiCGSTAB(     SparseMatrix<T> *A, T *b, T *x);
	int SolveBiCGSTAB2( SparseMatrix<T> *A, T *b, T *x);
	int SolveJacobiMethod( SparseMatrix<T> *A, T *b, T *x);
	int SolveGaussSeidelMethod( SparseMatrix<T> *A, T *b, T *x);
	int SolveSuccessiveOverRelaxationMethod( SparseMatrix<T> *A, T *b, T *x, T omega);
	int SolveExplicit( SparseMatrix<T> *A, T *b, T *x);
	int SolveLowerTriangular( SparseMatrix<T> *A, T *b, T *x);
	int SolveUpperTriangular( SparseMatrix<T> *A, T *b, T *x);
	int SolveLU( SparseMatrix<T> *LU, T *b, T *y, T *x);
	//int SolveBiCGSTAB( SparseMatrix *A, float *b, float *x);

	void PreconditionJacobi( SparseMatrix<T> *A, T *b);
	void PreconditionSOR( SparseMatrix<T> *A, T *b, T omega);
	void SuccessiveOverRelaxationPreconditioner( SparseMatrix<T> *A, SparseMatrix<T> *M);
	void IncompleteLUfactorization( SparseMatrix<T> *A);
	void IncompleteLUfactorization( SparseMatrix<T> *A, SparseMatrix<T> *M);
	//void PreconditionJacobi( SparseMatrix *A, float *b);
	//void PreconditionJacobi( SparseMatrix *A, double *b);
};

/*template <class T> void calculate_x( T *x, T omega, T *w, T alpha, T *p, int N);
template <class T> void calculate_r( T *r, T omega, T *v, T alpha, T *s, int N);
template <class T> void calculate_p( T *x, T *y, T alpha, T *a, T beta, T *b, int N);
template <class T> void calculate_p1( T *x, T *y, T alpha, T *a, int N);*/
void calculate_x( float *x, float omega, float *w, float alpha, float *p, int N);
void calculate_r( float *r, float omega, float *v, float alpha, float *s, int N);
void calculate_p( float *x, float *y, float alpha, float *a, float beta, float *b, int N);
void calculate_p1( float *x, float *y, float alpha, float *a, int N);
void calculate_x( double *x, double omega, double *w, double alpha, double *p, int N);
void calculate_r( double *r, double omega, double *v, double alpha, double *s, int N);
void calculate_p( double *x, double *y, double alpha, double *a, double beta, double *b, int N);
void calculate_p1( double *x, double *y, double alpha, double *a, int N);

#include "Solver.tcc"


#endif /* SOLVER_H_ */
